This script is used to infer trajectory from HF-SCs, using TInGa.
library(dplyr)
library(patchwork)
library(ggplot2)
.libPaths()
## [1] "/usr/local/lib/R/library"
In this section, we set the global settings of the analysis. We will store data there :
save_name = "non_matrix"
out_dir = "."
We load the sample information :
sample_info = readRDS(paste0(out_dir, "/../../1_metadata/hs_hd_sample_info.rds"))
project_names_oi = sample_info$project_name
graphics::pie(rep(1, nrow(sample_info)),
col = sample_info$color,
labels = sample_info$project_name)
This are the settings for trajectory inference :
traj_dimred = "harmony" # to infer trajectory on
traj_2d = "harmony_dm" # just of visualization
seed = 1337L
traj_max_dims = 40 # to infer trajectory on
name2D = "harmony_18_tsne" # just of visualization
We load the Seurat object :
sobj = readRDS(paste0(out_dir, "/", save_name, "_sobj.rds"))
sobj
## An object of class Seurat
## 17837 features across 5599 samples within 1 assay
## Active assay: RNA (17837 features, 2000 variable features)
## 8 dimensional reductions calculated: RNA_pca, RNA_pca_18_tsne, RNA_pca_18_umap, harmony, harmony_18_umap, harmony_18_tsne, harmony_dm, harmony_dm_5_umap
We set the root manually :
root_cell_id = "2022_14_GAACACTAGTGCCTCG-1"
sobj$is_root = colnames(sobj) == root_cell_id
sobj$cell_name = colnames(sobj)
root_plot = aquarius::plot_label_dimplot(sobj, reduction = traj_2d,
col_by = "is_root", col_color = c("gray80", "red"),
label_by = "cell_name", label_val = root_cell_id) +
ggplot2::ggtitle("Root cell") +
ggplot2::theme(plot.title = element_text(hjust = 0.5, face = "bold"),
aspect.ratio = 1) +
Seurat::NoLegend() + Seurat::NoAxes()
sobj$is_root = NULL
sobj$cell_name = NULL
root_plot
Now, we perform the trajectory inference :
set.seed(seed)
my_traj = aquarius::traj_tinga(sobj,
seed = seed,
expression_assay = "RNA",
expression_slot = "data",
count_assay = "RNA",
count_slot = "counts",
dimred_name = traj_dimred,
dimred_max_dim = traj_max_dims,
root_cell_id = root_cell_id,
tinga_parameters = list(max_nodes = 3))
## Add pseudotime
sobj$pseudotime = my_traj$pseudotime
(Time to run : 9.31 s)
We visualize pseudotime :
Seurat::FeaturePlot(sobj, features = "pseudotime",
reduction = traj_2d,
cols = viridis::viridis(n = 100)) +
ggplot2::lims(x = range(sobj@reductions[[traj_2d]]@cell.embeddings[, 1]),
y = range(sobj@reductions[[traj_2d]]@cell.embeddings[, 2])) +
ggplot2::labs(title = "Pseudotime") +
ggplot2::theme(aspect.ratio = 1,
plot.title = element_text(hjust = 0.5)) +
Seurat::NoAxes()
We visualize the trajectory with arrows :
dynplot::plot_dimred(trajectory = my_traj,
label_milestones = FALSE,
size_milestones = 0,
plot_trajectory = TRUE,
dimred = sobj[[traj_2d]]@cell.embeddings,
color_cells = 'pseudotime', color_trajectory = "none")
We color cells according to cluster ID :
dynplot::plot_dimred(trajectory = my_traj,
label_milestones = FALSE,
plot_trajectory = TRUE,
dimred = sobj[[traj_2d]]@cell.embeddings,
grouping = sobj$seurat_clusters,
size_milestones = 0,
color_trajectory = "none")
We visualize pseudotime :
Seurat::FeaturePlot(sobj, features = "pseudotime",
reduction = name2D,
cols = viridis::viridis(n = 100)) +
ggplot2::lims(x = range(sobj@reductions[[name2D]]@cell.embeddings[, 1]),
y = range(sobj@reductions[[name2D]]@cell.embeddings[, 2])) +
ggplot2::labs(title = "Pseudotime") +
ggplot2::theme(aspect.ratio = 1,
plot.title = element_text(hjust = 0.5)) +
Seurat::NoAxes()
We visualize the trajectory with arrows :
dynplot::plot_dimred(trajectory = my_traj,
label_milestones = FALSE,
size_milestones = 0,
plot_trajectory = TRUE,
dimred = sobj[[name2D]]@cell.embeddings,
color_cells = 'pseudotime', color_trajectory = "none")
We color cells according to cluster ID :
dynplot::plot_dimred(trajectory = my_traj,
label_milestones = FALSE,
plot_trajectory = TRUE,
dimred = sobj[[name2D]]@cell.embeddings,
grouping = sobj$seurat_clusters,
size_milestones = 0,
color_trajectory = "none")
We save the Seurat object :
saveRDS(sobj, file = paste0(out_dir, "/", save_name, "_sobj_traj_tinga.rds"))
We save the trajectory object :
saveRDS(my_traj, file = paste0(out_dir, "/", save_name, "_my_traj_tinga.rds"))
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.6 LTS
##
## Matrix products: default
## BLAS: /usr/local/lib/R/lib/libRblas.so
## LAPACK: /usr/local/lib/R/lib/libRlapack.so
##
## locale:
## [1] C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] dynutils_1.0.5 dynwrap_1.2.1 purrr_0.3.4 tidyr_1.1.4
## [5] ggplot2_3.3.5 patchwork_1.1.2 dplyr_1.0.7
##
## loaded via a namespace (and not attached):
## [1] softImpute_1.4 graphlayouts_0.7.0
## [3] pbapply_1.4-2 lattice_0.20-41
## [5] haven_2.3.1 dyndimred_1.0.3
## [7] vctrs_0.3.8 usethis_2.0.1
## [9] blob_1.2.1 survival_3.2-13
## [11] prodlim_2019.11.13 later_1.3.0
## [13] DBI_1.1.1 R.utils_2.11.0
## [15] SingleCellExperiment_1.8.0 rappdirs_0.3.3
## [17] uwot_0.1.8 dqrng_0.2.1
## [19] gng_0.1.0 jpeg_0.1-8.1
## [21] zlibbioc_1.32.0 pspline_1.0-18
## [23] pcaMethods_1.78.0 mvtnorm_1.1-1
## [25] htmlwidgets_1.5.4 GlobalOptions_0.1.2
## [27] future_1.22.1 UpSetR_1.4.0
## [29] laeken_0.5.2 leiden_0.3.3
## [31] clustree_0.4.3 lmds_0.1.0
## [33] parallel_3.6.3 scater_1.14.6
## [35] irlba_2.3.3 DEoptimR_1.0-9
## [37] tidygraph_1.1.2 Rcpp_1.0.9
## [39] readr_2.0.2 KernSmooth_2.23-17
## [41] carrier_0.1.0 promises_1.1.0
## [43] gdata_2.18.0 DelayedArray_0.12.3
## [45] limma_3.42.2 pkgload_1.2.2
## [47] graph_1.64.0 RcppParallel_5.1.4
## [49] Hmisc_4.4-0 fs_1.5.2
## [51] RSpectra_0.16-0 fastmatch_1.1-0
## [53] ranger_0.12.1 digest_0.6.25
## [55] png_0.1-7 sctransform_0.2.1
## [57] cowplot_1.0.0 DOSE_3.12.0
## [59] here_1.0.1 TInGa_0.0.0.9000
## [61] dynplot_1.1.0 ggraph_2.0.3
## [63] pkgconfig_2.0.3 GO.db_3.10.0
## [65] DelayedMatrixStats_1.8.0 gower_0.2.1
## [67] ggbeeswarm_0.6.0 iterators_1.0.12
## [69] DropletUtils_1.6.1 reticulate_1.26
## [71] clusterProfiler_3.14.3 SummarizedExperiment_1.16.1
## [73] circlize_0.4.15 beeswarm_0.4.0
## [75] GetoptLong_1.0.5 xfun_0.35
## [77] bslib_0.3.1 zoo_1.8-10
## [79] tidyselect_1.1.0 GA_3.2
## [81] reshape2_1.4.4 ica_1.0-2
## [83] pcaPP_1.9-73 viridisLite_0.3.0
## [85] rtracklayer_1.46.0 rlang_1.0.2
## [87] hexbin_1.28.1 jquerylib_0.1.4
## [89] dyneval_0.9.9 glue_1.4.2
## [91] waldo_0.3.1 RColorBrewer_1.1-2
## [93] matrixStats_0.56.0 stringr_1.4.0
## [95] lava_1.6.7 europepmc_0.3
## [97] DESeq2_1.26.0 recipes_0.1.17
## [99] labeling_0.3 httpuv_1.5.2
## [101] class_7.3-17 BiocNeighbors_1.4.2
## [103] DO.db_2.9 annotate_1.64.0
## [105] jsonlite_1.7.2 XVector_0.26.0
## [107] bit_4.0.4 mime_0.9
## [109] aquarius_0.1.5 Rsamtools_2.2.3
## [111] gridExtra_2.3 gplots_3.0.3
## [113] stringi_1.4.6 processx_3.5.2
## [115] gsl_2.1-6 bitops_1.0-6
## [117] cli_3.0.1 batchelor_1.2.4
## [119] RSQLite_2.2.0 randomForest_4.6-14
## [121] data.table_1.14.2 rstudioapi_0.13
## [123] org.Mm.eg.db_3.10.0 GenomicAlignments_1.22.1
## [125] nlme_3.1-147 qvalue_2.18.0
## [127] scran_1.14.6 locfit_1.5-9.4
## [129] scDblFinder_1.1.8 listenv_0.8.0
## [131] ggthemes_4.2.4 gridGraphics_0.5-0
## [133] R.oo_1.24.0 dbplyr_1.4.4
## [135] BiocGenerics_0.32.0 TTR_0.24.2
## [137] readxl_1.3.1 lifecycle_1.0.1
## [139] timeDate_3043.102 ggpattern_0.3.1
## [141] munsell_0.5.0 cellranger_1.1.0
## [143] R.methodsS3_1.8.1 proxyC_0.1.5
## [145] visNetwork_2.0.9 caTools_1.18.0
## [147] codetools_0.2-16 Biobase_2.46.0
## [149] GenomeInfoDb_1.22.1 vipor_0.4.5
## [151] lmtest_0.9-38 msigdbr_7.5.1
## [153] htmlTable_1.13.3 triebeard_0.3.0
## [155] lsei_1.2-0 xtable_1.8-4
## [157] ROCR_1.0-7 BiocManager_1.30.10
## [159] scatterplot3d_0.3-41 abind_1.4-5
## [161] farver_2.0.3 parallelly_1.28.1
## [163] RANN_2.6.1 askpass_1.1
## [165] GenomicRanges_1.38.0 RcppAnnoy_0.0.16
## [167] tibble_3.1.5 ggdendro_0.1-20
## [169] cluster_2.1.0 future.apply_1.5.0
## [171] Seurat_3.1.5 dendextend_1.15.1
## [173] Matrix_1.3-2 ellipsis_0.3.2
## [175] prettyunits_1.1.1 lubridate_1.7.9
## [177] ggridges_0.5.2 igraph_1.2.5
## [179] RcppEigen_0.3.3.7.0 fgsea_1.12.0
## [181] remotes_2.4.2 scBFA_1.0.0
## [183] destiny_3.0.1 VIM_6.1.1
## [185] testthat_3.1.0 htmltools_0.5.2
## [187] BiocFileCache_1.10.2 yaml_2.2.1
## [189] utf8_1.1.4 plotly_4.9.2.1
## [191] XML_3.99-0.3 ModelMetrics_1.2.2.2
## [193] e1071_1.7-3 foreign_0.8-76
## [195] withr_2.5.0 fitdistrplus_1.0-14
## [197] BiocParallel_1.20.1 xgboost_1.4.1.1
## [199] bit64_4.0.5 foreach_1.5.0
## [201] robustbase_0.93-9 Biostrings_2.54.0
## [203] GOSemSim_2.13.1 data.tree_1.0.0
## [205] rsvd_1.0.3 memoise_2.0.0
## [207] evaluate_0.18 forcats_0.5.0
## [209] rio_0.5.16 geneplotter_1.64.0
## [211] tzdb_0.1.2 caret_6.0-86
## [213] ps_1.6.0 DiagrammeR_1.0.6.1
## [215] curl_4.3 fdrtool_1.2.15
## [217] fansi_0.4.1 highr_0.8
## [219] urltools_1.7.3 xts_0.12.1
## [221] GSEABase_1.48.0 acepack_1.4.1
## [223] edgeR_3.28.1 checkmate_2.0.0
## [225] scds_1.2.0 cachem_1.0.6
## [227] randomForestSRC_2.12.1 desc_1.4.1
## [229] npsurv_0.4-0 babelgene_22.3
## [231] rjson_0.2.20 openxlsx_4.1.5
## [233] ggrepel_0.9.1 clue_0.3-60
## [235] rprojroot_2.0.2 stabledist_0.7-1
## [237] tools_3.6.3 sass_0.4.0
## [239] nichenetr_1.1.1 magrittr_2.0.1
## [241] RCurl_1.98-1.2 proxy_0.4-24
## [243] car_3.0-11 ape_5.3
## [245] ggplotify_0.0.5 xml2_1.3.2
## [247] httr_1.4.2 assertthat_0.2.1
## [249] rmarkdown_2.18 boot_1.3-25
## [251] globals_0.14.0 R6_2.4.1
## [253] Rhdf5lib_1.8.0 nnet_7.3-14
## [255] RcppHNSW_0.2.0 progress_1.2.2
## [257] genefilter_1.68.0 statmod_1.4.34
## [259] gtools_3.8.2 shape_1.4.6
## [261] HDF5Array_1.14.4 BiocSingular_1.2.2
## [263] rhdf5_2.30.1 splines_3.6.3
## [265] AUCell_1.8.0 carData_3.0-4
## [267] colorspace_1.4-1 generics_0.1.0
## [269] stats4_3.6.3 base64enc_0.1-3
## [271] dynfeature_1.0.0 smoother_1.1
## [273] gridtext_0.1.1 pillar_1.6.3
## [275] tweenr_1.0.1 sp_1.4-1
## [277] ggplot.multistats_1.0.0 rvcheck_0.1.8
## [279] GenomeInfoDbData_1.2.2 plyr_1.8.6
## [281] gtable_0.3.0 zip_2.2.0
## [283] knitr_1.41 ComplexHeatmap_2.14.0
## [285] latticeExtra_0.6-29 biomaRt_2.42.1
## [287] IRanges_2.20.2 fastmap_1.1.0
## [289] ADGofTest_0.3 copula_1.0-0
## [291] doParallel_1.0.15 AnnotationDbi_1.48.0
## [293] vcd_1.4-8 babelwhale_1.0.1
## [295] openssl_1.4.1 scales_1.1.1
## [297] backports_1.2.1 S4Vectors_0.24.4
## [299] ipred_0.9-12 enrichplot_1.6.1
## [301] hms_1.1.1 ggforce_0.3.1
## [303] Rtsne_0.15 shiny_1.7.1
## [305] numDeriv_2016.8-1.1 polyclip_1.10-0
## [307] grid_3.6.3 lazyeval_0.2.2
## [309] Formula_1.2-3 tsne_0.1-3
## [311] crayon_1.3.4 MASS_7.3-54
## [313] pROC_1.16.2 viridis_0.5.1
## [315] dynparam_1.0.0 rpart_4.1-15
## [317] zinbwave_1.8.0 compiler_3.6.3
## [319] ggtext_0.1.0